
from libatomism import *
from math import *
import array 
import matplotlib.pyplot as plt

msLogger.setPriorityMin(Priority.ERROR)

unit = msUnitsManager.New("Angstrom Degree amu kcal/mol ps")


zmat = msZMat.New(unit).setId("zmat").set("C C 0 1.52 O 1 1.43 0 108.0 H 0 1.09 1 111 2 180 H 0 1.09 3 109 1 121. H 0 1.09 4 108 3 -118 H 2 0.97 1 108 0 180 H 1 1.1 2 111 0 120 H 1 1.1 2 111 0 -120 ")

system = msSystem.New(unit).addEntity(zmat)

gcoors = msGeneralizedCoordinates.New(unit)
q0     = msScalarVariable.New("Angstrom").set(0., 0.,  180., 0.01, 0.1).setId("q0")
q1     = msScalarVariable.New("Angstrom").set(0., 0.,  60.,  0.01, 0.1).setId("q1")

gcoors.addVar(q0).addVar(q1)

kinop = msKineticOperator.New(unit).set(gcoors,system)
kinop.addDynamicDofLinear( zmat, 14, q0).addDynamicDofLinear( zmat, 5, q1)
print kinop
viewer= msVisualizerAtoms.New()
viewer.watch(system)

cfour = qcCfour.New().setLevel( BasisSet.ANO0 , CiType.CCSDpT ).setGenBasSymLink("/org/centers/pecos/APPLICATIONS/cfour/basis/GENBAS")

emin    = cfour.singlePoint(zmat)
print "emin=",emin

file = open('dataAI.txt','w')

for i in range(100):
    gcoors.setRandom()
    kinop.setDynamicDofs()
    e = cfour.singlePoint(zmat)
    file.write(str(q0.Value)+" "+str(q1.Value)+" "+str(e)+" "+str(unit.convert("hartree",e-emin))+"\n" )
    print q0.value," ",q1.Value," ",e,unit.convert("hartree",e-emin)
file.close()
  
